/*
fucntion: non-maximum suppression
input:
pMag:   pointer to Magnitude,
pGradX: gradient of x-direction
pGradY: gradient of y-direction
sz: size of pMag (width = size.cx, height = size.cy)
output:
pNSRst: result of non-maximum suppression

https://blog.csdn.net/kezunhai/article/details/11620357
*/
void NonMaxSuppress(int *pMag, int *pGradX, int *pGradY, size_t sz, unsigned char *pNSRst)
{
    long x, y;
    int nPos;
    // the component of the gradient
    int gx, gy;
    // the temp varialbe
    int g1, g2, g3, g4;
    double weight;
    double dTemp, dTemp1, dTemp2;
    // 设置图像边缘为不可能的分界点
    for (x = 0; x < sz.cx; x++)
    {
        pNSRst[x] = 0;
        pNSRst[(sz.cy - 1) * sz.cx + x] = 0;
    }
    for (y = 0; y < sz.cy; y++)
    {
        pNSRst[y * sz.cx] = 0;
        pNSRst[y * sz.cx + sz.cx - 1] = 0;
    }

    for (y = 1; y < sz.cy - 1; y++)
    {
        for (x = 1; x < sz.cx - 1; x++)
        {
            nPos = y * sz.cx + x;
            // if pMag[nPos]==0, then nPos is not the edge point
            if (pMag[nPos] == 0)
            {
                pNSRst[nPos] = 0;
            }
            else
            {
                // the gradient of current point
                dTemp = pMag[nPos];
                // x,y 方向导数
                gx = pGradX[nPos];
                gy = pGradY[nPos];
                // 如果方向导数y分量比x分量大，说明导数方向趋向于y分量
                if (abs(gy) > abs(gx))
                {
                    // calculate the factor of interplation
                    weight = fabs(gx) / fabs(gy);
                    g2 = pMag[nPos - sz.cx]; // 上一行
                    g4 = pMag[nPos + sz.cx]; // 下一行
                    // 如果x,y两个方向导数的符号相同
                    // C 为当前像素，与g1-g4 的位置关系为：
                    // g1 g2
                    //    C
                    //    g4 g3
                    if (gx * gy > 0)
                    {
                        g1 = pMag[nPos - sz.cx - 1];
                        g3 = pMag[nPos + sz.cx + 1];
                    }
                    // 如果x,y两个方向的方向导数方向相反
                    // C是当前像素，与g1-g4的关系为：
                    //     g2 g1
                    //     C
                    //  g3 g4
                    else
                    {
                        g1 = pMag[nPos - sz.cx + 1];
                        g3 = pMag[nPos + sz.cx - 1];
                    }
                }
                else
                {
                    // 插值比例
                    weight = fabs(gy) / fabs(gx);
                    g2 = pMag[nPos + 1]; // 后一列
                    g4 = pMag[nPos - 1]; // 前一列
                                         // 如果x,y两个方向的方向导数符号相同
                    // 当前像素C与 g1-g4的关系为
                    //  g3
                    //  g4 C g2
                    //        g1
                    if (gx * gy > 0)
                    {
                        g1 = pMag[nPos + sz.cx + 1];
                        g3 = pMag[nPos - sz.cx - 1];
                    }

                    // 如果x,y两个方向导数的方向相反
                    //  C与g1-g4的关系为
                    //  g1
                    //  g4 C g2
                    //       g3
                    else
                    {
                        g1 = pMag[nPos - sz.cx + 1];
                        g3 = pMag[nPos + sz.cx - 1];
                    }
                }

                dTemp1 = weight * g1 + (1 - weight) * g2;
                dTemp2 = weight * g3 + (1 - weight) * g4;
                // 当前像素的梯度是局部的最大值
                // 该点可能是边界点
                if (dTemp >= dTemp1 && dTemp >= dTemp2)
                {
                    pNSRst[nPos] = 128;
                }
                else
                {
                    // 不可能是边界点
                    pNSRst[nPos] = 0;
                }
            }
        }
    }
}
